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Abstract 

We report ab initio calculations of the melting curve of molybdenum for the pressure range 
— 400 GPa. The calculations employ density functional theory (DFT) with the Perdew-Burke- 
Ernzerhof exchange-correlation functional in the projector augmented wave (PAW) implementa- 
tion. We present tests showing that these techniques accurately reproduce experimental data on 
low-temperature b.c.c. Mo, and that PAW agrees closely with results from the full-potential lin- 
earized augmented plane-wave implementation. The work attempts to overcome the uncertainties 
inherent in earlier DFT calculations of the melting curve of Mo, by using the "reference coexistence" 
technique to determine the melting curve. In this technique, an empirical reference model (here, 
the embedded-atom model) is accurately fitted to DFT molecular dynamics data on the liquid and 
the high-temperature solid, the melting curve of the reference model is determined by simulations 
of coexisting solid and liquid, and the ab initio melting curve is obtained by applying free-energy 
corrections. Our calculated melting curve agrees well with experiment at ambient pressure and is 
consistent with shock data at high pressure, but does not agree with the high pressure melting 
curve deduced from static compression experiments. Calculated results for the radial distribution 
function show that the short-range atomic order of the liquid is very similar to that of the high-T 
solid, with a slight decrease of coordination number on passing from solid to liquid. The electronic 
densities of states in the two phases show only small differences. The results do not support a 
recent theory according to which very low dT m /dP values are expected for b.c.c. transition metals 
because of electron redistribution between s-p and d states. 
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I. INTRODUCTION 



Over the past five years, a major controversy has developed about the high-pressure 
melting curves of transition metals. According to static compression experiments, performed 
using diamond anvil cells (DAC) at pressures from ambient up to ~ 100 GPa (1 Mbar), the 
melting temperatures T m of many transition metals, particularly those having the b.c.c. 
crystal structure, change by no more than a few hundred K over this pressure ranged. In 
some cases, this finding seems to be in gross conflict with shock experiments, which indicate 
an increase of T m of several thousand K over the same pressure range^^. Theoretical 
work->£&2 based directly or indirectly on density functional theory (DFT) generally supports 
the shock data. We report here a detailed DFT study of the high-pressure melting curve of 
Mo, a metal for which there are some of the largest differences between DAC measurements 
and other data (see e.g. Fig. 2 in Ref.— ). 

There has already been quite extensive theoretical work on the high-P melting of metals. 
Some of this has been motivated by the desire to understand the properties of solid and 
liquid Fe in the Earth's coreii. DFT-based calculations of the Fe melting curve^^ 1 ^ 1 ^ 
up the pressure at the boundary between the solid inner core and the liquid outer core 
provide one of the important ways of constraining the temperature distribution in the core. 
Disagreements with DAC measurements on Fe and other transition metals are cause for 
concern, because if DFT were shown to be seriously in error, the reliability of DFT for the 
study of planetary interiors would be called into question. But even without this practical 
motivation, a major disagreement between DFT predictions and experimental data must 
be taken seriously, because it suggests an unexpected failure either of commonly used DFT 
approximations, or of apparently well established experimental techniques. Reassuringly, 
DFT melting curves agree very closely with experiment for some metals, including Al^^ 
and Cu^. The large disagreements arise mainly for transition metals, and the suggestion is 
that they are linked to d-band bonding. 

Since the melting slope dT m /dP is equal by the Clausius-Clapeyron relation to V ls /S ls , 
where V ls and S ls are the volume and entropy of fusion, and since S ls is unlikely to have 
exceptionally large values, a very low dT m /dP is likely to be due to a low volume of fusion. 
It has been arguec-i^ that a low V ls might be expected for b.c.c. transition metals, because 
the liquid may be more close packed than the solid. But recently, an additional argument 
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has been advanced for low dT m /dP values^. This argument depends on the fact that the 
distribution of conduction electrons between s-p states and d states is known to depend 
on the degree of compression and on the crystal stucture^ 1 ^. The suggestion is that the 
(assumed) change of coordination on going from b.c.c. solid to liquid leads to a change 
of electronic structure, and hence a change in the electron distribution between s-p and d 
states, and that this redistribution stabilises the liquid and lowers T m . One of the important 
purposes of the present work is to use DFT molecular dynamics (m.d.) simulations to test 
these suggestions for the case of Mo. 

There has been previous DFT work- 1 ^ on the high-pressure melting of Mo. The early 
work of Moriarty 6 employed a many-body total energy function derived from first-principles 
"generalized pseudopotential theory"—. This approximates the total energy of the system 
in terms of volume dependent 1-, 2-, 3- and 4-body interatomic potentials, and accounts 
for the angular forces that are known to be important in transition metals. This form of 
total-energy function is designed to be fully transferable between different structures, and 
should be valid for both the solid and the liquid state. The model total-energy function was 
used in m.d. simulations to determine the melting curve by two methods. The first method 
consisted of cycling the simulated system up and down through the melting point at fixed 
volume. The results were cross-checked against a second method based on the calculation of 
the free energies of the solid and the liquid. This was pioneering work, but its quantitative 
accuracy can be questioned, because some of the phonon frequencies predicted by the total 
energy function agreed rather poorly with experiment. In addition, the surprising claim 
was made that T m could be changed by up to a factor of two by the inclusion of thermal 
electronic excitations, which were treated only crudely. The predicted melting curve was 
consistent with shock data, but was far above the curve given by recent DAC measurements: 
the difference of T m values amounts to ~ 3000 K at P = 100 GPa. Much more recently, 
direct DFT m.d. simulation has been used by Belonoshko et air- to map out the Mo melting 
curve up to P ~ 330 GPa. Their melting curve, like that of Moriarty^, is consistent with 
shock data, but disagrees strongly with DAC data. However, even though the simulations 
employed an implementation of DFT that is expected to be accurate, the method used 
to map the melting curve may still give inaccurate results, as noted by the authors. Their 
method was to heat the simulated solid at constant volume until the internal energy, pressure, 
radial distribution function and self-diffusion coefficient showed discontinuities attributable 
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to melting. This approach gives an upper bound to T m , but may overestimate it significantly, 
because of superheating. The authors estimated that the error in T m due to superheating 
should be ~ 20 %, but the arguments for this estimate are indirect. 

The present new work on Mo melting has several aims. The first aim is to calculate 
more accurately the melting curve that follows from the adopted DFT exchange-correlation 
functional E xc . (We present tests showing that the GGA-PBE functional is a good choice.) 
It is only by doing this that the possible reasons for the disagreement between theory and 
experiment can be narrowed down. We have shown in our recent work on other materials 
that the errors in computed DFT T m values can be reduced to a few percent, and we have 
described several techniques for doing thisi 1 -. The technique used here is the "reference 
coexistence method" 24 . This requires the accurate fitting of an empirical "reference" total 
energy function to DFT m.d. simulations on the solid and the liquid; the melting curve of the 
reference model is then calculated from simulations on large systems consisting of coexisting 
solid and liquid; as an essential last step, free-energy corrections for the difference between 
the reference and DFT total-energy functions are used to correct the melting curve. This 
technique was successfully used in our work on the melting of Cu 1 ^, and it has been shown to 
give results in excellent agreement with an alternative technique, in which DFT free energies 
of the solid and liquid are calculated 1 ^ 2 ^. Since it is a 'thermodynamic' technique relying 
on equality of Gibbs free energies of the two phases, it cannot suffer from superheating 
problems. A second important aim is to study the differences of atomic and electronic 
structure of the coexisting solid and liquid, to provide an improved understanding of the 
factors that determine the melting curve of Mo. The empirical reference model used to 
determine the DFT melting curve has the form of the embedded atom model (EAM)- 1 ^. 
A useful side benefit of the work is that we obtained a parameterized EAM that mimics 
quite well the DFT total-energy function of high-T solid and liquid Mo. This model reveals 
important features of the energetics of Mo at high-P and high-T, and we expect it to be 
useful in future modelling work on this metal. 

The remainder of the paper is organized as follows. A brief summary of the DFT mod- 
elling techniques is given in Sec. [Til followed by an outline of the reference coexistence 
technique. Then, Sec. II III details the rather extensive tests we have performed to ensure 
that the techniques deliver an accurate description of the energetics and the vibrational and 
electronic properties of Mo. Sec. IIVI presents our results for the DFT melting curve of Mo 
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up to 400 GPa, together with the entropy and volume of fusion as a function of pressure; 
our comparison of the atomic and electronic structures of the high-T b.c.c. solid and the 
liquid is reported at the end of the Section. A discussion of all the results and their relation 
with previous work is given in Sec. |V] followed by our conclusions. 

II. TECHNIQUES 
A. DFT methods 

A comprehensive description of DFT methods as applied to the modelling of condensed 
matter is given in recent books^^i. As is well known, there is one and only one uncontrol- 
lable approximation in DFT, namely the approximation used for the exchange-correlation 
functional E xc . There is abundant evidence that commonly used E xc functionals yield accu- 
rate results for a range of properties of transition-metal crystals, including equilibrium lattice 
parameter, elastic constants, phonon frequencies, T = equation of state, solid-state phase 
boundaries, etc^i. DFT provides the basis for our present understanding of the electronic 
structure of transition-metal crystals as a function of pressure. There is also some evidence 
that it accurately reproduces the structure of liquid transition metals 2 '. Nevertheless, we 
considered it necessary to test the accuracy of different E xc in the case of Mo, as we report 
in Sec. EL 

A completely separate issue from the choice of E xc is the implementation of DFT that 
is used. This concerns mainly the way that the electron orbitals are represented. For 
simulations of the high-T solid and the liquid, DFT molecular dynamics (m.d.) must be 
used, and for this we use the PAW (projector augmented wave) technique 2 ^ 1 ^, which is 
generally regarded as the most accurate for m.d. purposes. The PAW method developed 
by Blochl 28 efficiently combines some of the features originally devised within both linear 
augmented-plane-wave (LAPW, briefly described below) and pseudopotential approaches. 
In PAW, a linear transformation between the all-electron and pseudized wavefunctions is 
defined in terms of all-electron and pseudized partial waves and a set of projector functions 
localized on the atoms. There is a simple formal relationship between PAW and the ultra- 
soft pseudopotential method of Vanderbilt^, and computationally the two approaches are 
almost equivalent^. 
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For m.d., we use the method known as Born-Oppenheimer dynamics, in which the self- 
consistent ground state is recalculated at each new m.d. time-step. Given the criticisms 
that are sometimes made of DFT work on melting properties^, is important to be clear 
that this means the recalculation of the electronic structure of the entire simulated system, 
so that changes of electronic structure on going from solid to liquid at any degree of com- 
pression are fully included. Although PAW is well suited to m.d. simulation, we regard it as 
essential to check its accuracy against still more accurate implementations of DFT, and to 
do this we have compared with the predictions of the FP-LAPW (full-potential linearized 
augmented plane wave) technique, which was developed for high-precision calculations on 
crystals^i^ 1 ^. In this approach, space is divided into spherical regions centred on the ions, 
and the interstitial space between the spheres. Accurate solution of Schrodinger's equation 
on a radial grid for each angular momentum within each sphere allows both valence and core 
states to be treated to arbitrarily high precision. The FP-LAPW results to be presented 
include fully relativistic treatment of the core states and scalar relativistic treatment of the 
states in valence^. In FP-LAPW calculations on crystalline Mo, we can ensure that they 
are fully converged with respect to all technical parameters, so that all errors of implemen- 
tation are negligible, and the only approximation is that due to E xc itself. All our PAW and 
FP-LAPW calculations were performed using the VASP code^ and the WIEN2k code^, 
respectively. 

Normally, the aim of DFT calculations is to obtain the electronic ground state for given 
ionic positions. However, because of the high T involved, it is essential in the present 
work to include thermal electronic excitations, so that for any given ionic positions we must 
determine the orbitals and occupation numbers that self-consistently minimize the electronic 
free energy, using the T > version of DFT originally developed by Mermin^. We know 
this is essential, because work on other transition metals™ shows that the electronic specific 
heat becomes comparable with the vibrational contribution for T ~ 5000 K. Without this, 
the free energy difference between solid and liquid might be seriously in error. All our m.d. 
simulations are done in the canonical (N, V, T) ensemble, with the electronic T set equal to 
the T of the ensemble. 

Values of technical parameters (fc-point sampling, plane-wave cut-off, etc.) will be given 
when we present the calculations. 
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B. Reference coexistence methods 



There are three steps in the reference coexistence technique^. First, an empirical refer- 
ence model is fitted to ab initio simulations of the solid and the liquid at thermodynamic 
states close to the expected melting curve. Then, the reference model is used to perform 
simulations on large systems in which solid and liquid coexist, so as to obtain points on 
the melting curve of the model. Finally, the differences between the reference and ab initio 
total energy functions are used to correct the melting properties of the fitted model so as to 
obtain the ab initio melting properties. 

In this work, we have used the embedded-atom model (EAM)^^ as the reference model. 
The total energy function U Te [ of this model for a system of iV atoms has the form: 



where = | — r 3 - 1 is the separation of atoms i and j. The first term on the right represents 
an inverse-power repulsive pair potential, while the second (embedding) term describes the 
(i-band bonding. The model is specified by the characteristic length a, the energy scale e, 
the dimensionless coefficient C characterizing the strength of the embedding energy, and the 
embedding and repulsive exponents m and n. We wish [7 re f (r 1; . . . rjy) to mimic as closely as 
possible the ab initio total energy function Uai(ti, ...!•#), and we achieve this by adjusting 
the EAM parameters. 

The method for fitting U re f to Uai is designed so as to minimize the corrections to the 
melting curve caused by the difference AU = Uai — U re {. We therefore recall the correction 
scheme before describing the fitting itself. For given P and T, the difference G^i = G l Al — 
G s Al between the Gibbs free energies of the ab initio liquid and solid deviates from the 
corresponding difference G\ s c{ = G Tei — G* ef of the reference liquid and solid, and we write: 



It is the shift AG ls (P, T) caused by changing the total-energy function from C/ re f to Uai that 
causes the shift of melting temperature at given pressure. To first order, the latter shift is^: 




m -I 1/2 



(1) 



G 1 a\{P, T) = G l T s c{ (P, T) + AG ls (P, T) . 



(2) 



V 

m 



A Qls ( T ref) 



(3) 
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where S^ f is the difference between the entropies of liquid and solid (i.e. the entropy of 
fusion) of the reference system, and AG ls (T™ ) is AG ls evaluated at the melting temperature 
of the reference system. 

The shift AG ls is the difference of shifts of Gibbs free energies of liquid and solid caused 
by the shift AU = Uai — U re f of total energy function. We find it convenient to perform 
our simulations at constant volume and temperature. Under these conditions, the shift of 
Helmholtz free energy AF arising from AU is given by the well-known expansion: 

AF = (AU) Iei - ^(3(5AU 2 ) rcf + ■■■ , (4) 

where (3 = l/k^T, 5AU = AU — (A£/) re f, and the averages are taken in the reference 
ensemble. From AF, we obtain the shift of Gibbs free energy at constant pressure as: 

AG = AF - ^Vk t (AP) 2 , (5) 

where kt is the isothermal compressibility and AP is the change of pressure when U re { is 
replaced by Uai at constant V and T. 

In fitting the EAM U re { to the ab initio Uai so as to minimize the corrections we have 
just described, we see from Eq. (jl]) and (jSj) that the effects of (AC/) re f, (5AU 2 ) re { and AP all 
need to be made small. Concerning (AU) re {, we note that addition of a position-independent 
constant to either Uai or U Te f has no effect on the properties of the solid or the liquid, since 
it simply redefines the energy zero, so that a large value of (AC/) re f is not in itself significant. 
However, it is crucially important that (AC/) re f should have almost the same value in the solid 
and the liquid, because otherwise there would be a large shift of AG ls . We also seek to make 
(5AU 2 ) re { and AP small in both phases. To satisfy all these requirements, we perform long ab 
initio m.d. simulations of the solid and the liquid at thermodynamic states near the expected 
melting curve. A large number of statistically independent configurations are then drawn 
from both these runs, so that we get a large set of Uai values for a collection of configurations 
representative of the solid and the liquid. Our procedure is then to minimize the mean square 
fluctuations of AU over this whole set: with AU the mean value of AU over the whole set, 
we minimize d~AU 2 , where SAU = AU — AU. Note that this 5AU 2 is not the same quantity 
as the {5AU 2 ) rc { appearing in Eq. (TJJ, because 8AU 2 characterizes the fluctuations of AU 
over a set of configurations drawn from both liquid and solid. Minimization of 5AU 2 has 
the effect of reducing simultaneously the difference of (AC/) re f between solid and liquid, and 
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of reducing (<5A£7 2 ) re f in the two phases. In order to reduce also AP 2 in the two phases, we 
add AP 2 into the quantity to be minimized, with a suitable weight. As will be described in 
Sec. IIV A[ we find it necessary to refit the EAM in this way for different ranges of pressure 
along the melting curve. A convenient way of characterizing the quality of fit of U rc { to 



Uai is obtained by dividing 5AU 2 by the mean square value of the fluctuation of ab initio 



energy 8U ai = Uai — Uai, where Uai is the ab initio energy averaged over the collection 
of 



of solid and liquid configurations. The dimensionless quantity characterizing the fit is then 

1/2 



5AU 2 /5U\ l ; the smaller the value of ?/>, the better is the fit. 
Following our previous work— our simulations of coexisting solid and liquid with 
the reference model are performed in the (N, V, E) ensemble. Starting with a supercell 
containing the perfect b.c.c. crystal, we thermalize it at a temperature slightly below the 
expected melting temperature. The system remains in the solid state. The simulation is then 
halted, and the positions of the atoms in one half of the supercell are held fixed, while the 
other half is heated to a very high temperature (typically eight times the expected melting 
temperature), so that it melts completely. With the fixed atoms still fixed, the molten part 
is then rethermalized to the expected melting temperature. Finally, the fixed atoms are 
released, thermal velocities are assigned, and the whole system is allowed to evolve freely 
at constant (N,V,E) for a long time (normally more than 60 ps), so that the solid and 
liquid come into equilibrium. The system is monitored by calculating the average number 
of particles in slices of the cell taken parallel to the boundary between the solid and liquid. 
With this protocol, there is a certain amount of trial and error to find the overall volume that 
yields the coexisting solid and liquid system. (An example of the density profile obtained 
when the two phases are in stable coexistence will be given in Sec. II VI ) Our main simulations 
were done on cells containing 6750 atoms, constructed as a 15 x 15 x 30 b.c.c. supercell, the 
long axis being perpendicular to the initial liquid-solid boundary. We tested the adequacy 
of this system size by repeating the simulations with larger systems (around 10000 atoms), 
and found no change in the results. 

Finally, our best value of the ab initio melting temperature is obtained by adding to the 
reference melting temperature the correction given by Eq. ([3]). The reference entropy 
of fusion S l r l f needed in this equation is obtained by performing independent reference m.d. 
simulations of the solid and the liquid in the (N, V, T) ensemble at the temperatures at 
which we made the coexistence simulations, using at each of these temperatures the solid 
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and liquid volumes that yield the coexistence pressure. These simulations give the enthalpy 
of fusion i^ref' fr° m which we obtain S l r l f from the relation H\ s ei = T^fSj^. To obtain the 
values of AG ls (T^) needed in Eq. © , we again perform reference m.d. simulations of solid 
and liquid at the coexistence temperatures and volumes, and calculate Uai for a statistically 
independent set of configurations drawn from these simulations. These Uai values must, of 
course, be fully converged with respect to system size and k-point sampling. We find that 
for a system of 125 atoms and with Uai calculated with a 2 x 2 x 2 k-point grid for both 
solid and liquid, Uai is converged within 5 meV/atom. These Uai values are then used to 
compute (AZ7)ref and (5AU 2 ) TC f . The same calculations yield AP for solid and liquid, which 
we use in Eq. (j5J). 

III. TESTS OF METHODS: ZERO-TEMPERATURE 

To assess the accuracy of the techniques, we have carried out a number of tests on the zero- 
temperature crystal in the pressure range — 300 GPa. There are three important questions 
to analyze. First, we want to test the accuracy of different exchange-correlation functionals 
compared with experimental data. Second, we aim to study the effect on PAW results of 
including different electronic states in the valence set. Third, FP-LAPW calculations are 
performed to assess errors incurred by the PAW approximation. In addition, we have done 
tests to demonstrate that PAW correctly reproduces the changes of electronic structure with 
pressure given by the FP-LAPW method. Tests of phonon frequencies are also carried out, 
because of their relevance to the vibrational free energy of the system. 

A. PAW and FP-LAPW calculations 

At very high pressures, states that would normally be treated as core states may respond 
significantly to compression. In the case of Mo, the 4p states lie only ~ 35 eV below 
the Fermi energy Ep and we always include them in the valence set. The 4s states lie 
considerably deeper at ~ 61 eV below E-p, and we have examined the effect of including 
them. We show in Fig. [Tj PAW results for the pressure as a function of volume at T = K, 
with and without the 4s states included in the valence set, compared with experimental 
data 5 . Results are shown with both LDA (Ceperley-Alder parameterization 37 ) and GGA 
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(Perdew-Burke-Ernzerhof form^) exchange-correlation functional. All calculations were 
performed on a primitive b.c.c. cell, with a 32 x 32 x 32 k-point grid and with energy 
cut-offs of 224.5 eV (without 4s states) and 287.6 eV (with 4s states); these settings ensure 
energy convergence to better than 1 meV/atom. With LDA, inclusion of 4s states makes 
a significant difference, and improves the agreement with experiment at high pressure, but 
with GGA the effect of including 4s states is very small. All the approximations deviate 
noticeably from experiment. At low pressures, LDA underestimates the volume by about 
3.1 %, while GGA overestimates it by 1.2 %. At high pressures, P ~ 300 GPa, LDA volumes 
are only 0.6 % below measured values, while GGA continues to overestimate them by about 
1.7 %. 

To check the accuracy of PAW itself, we have compared with FP-LAPW calculations of 
the T = K P(V) relation, performed with the WIEN2k codeB. Our FP-LAPW calculations 
include fully relativistic effects of the core states and scalar relativistic treatment of the states 
in valence^, and the tolerances are chosen so as to ensure energy convergence to better than 
1 meV/atom^. As shown in Fig. [2} the FP-LAPW P(V) results are almost indistinguishable 
from the PAW results over the entire pressure range of interest, with both LDA and GGA 
functionals. This provides valuable confirmation of the accuracy of the PAW technique, on 
which all our present calculations of melting properties are based. Out of interest, we also 
performed FP-LAPW calculations using the recently developed Wu-Cohen form of GGA—, 
which has been reported to give improved predictions of condensed-matter properties. This 
functional satisfies the same constraints used to construct the PBE functional, but the 
enhancement factor entering the exchange energy density is chosen to match closely the 
gradient expansion of Svendsen and von Barth^ for systems having a slowly varying density. 
For this reason, this functional is expected to perform well for solids but not so well for atoms 
or molecules, where the variation of the valence electron density is generally larger. We find 
that this functional is in much better agreement with the experimental P(V) data than 
either LDA or GGA(PBE). Unfortunately, at the time most of the present work was done, 
we had no way of performing Wu-Cohen calculations within PAW, so we were unable to do 
melting calculations with this new functional. 
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B. Electronic density of states 



The electronic structure of transition metals has been intensively studied, and changes 
of electronic structure with increasing pressure have been thoroughly investigated both ex- 
perimentally and theoretically^^ 1 ^ 1 ^. The energetics of these metals is, of course, strongly 
dominated by the d-bands, but the s-p bands also play an important role. An effect that 
has been well recognised for many years is the pressure induced upward shift of the s-p band 
relative to the d band, caused by the greater spatial extent of the s-p orbitals 21 . This relative 
shift means that increasing pressure causes a transfer of electrons from s-p to d bands, re- 
sulting in an increase of occupancy of the latter. It has been proposed recently^ 1 ^ that the 
large differences between theoretical predictions of melting curves and the results of static 
compression experiments for some transition metals may be due to an incorrect treatment 
of s-p — > d transfer in the DFT-based simulations. In order to demonstrate that our PAW 
calculations reproduce the electronic structure as given by the most accurate DFT methods, 
we have made detailed comparisons of the electronic density of states (DOS) calculated with 
PAW and FP-LAPW over the pressure range — 300 GPa. 

We show in Fig. [3] the electronic DOS from FP-LAPW and PAW calculated at 0, 150 
and 300 GPa, using the GGA(PBE) functional. We note the essentially perfect agreement 
between the two methods. At all pressures, the DOS consists of two separate parts: a narrow 
peak at about 35 eV below the Fermi energy, corresponding to 4p states, and a much broader 
distribution consisting of the multiple peaks due to the 4d bands superimposed on the slowly 
varying DOS of the s-p bands. We note the characteristic feature of b.c.c. transition metals 
that the Fermi energy falls in a deep minimum in the d-band DOS. As expected, the widths 
of the 4p and 4d parts of the DOS broaden markedly with increasing pressure. The pressure 
induced relative shift of s-p and d bands cannot be clearly seen from the DOS itself. However, 
it is very clear from the band structure, shown at pressures of and 300 GPa in Fig. HI The 
state lying ~ 7 eV below the Fermi energy is the bottom of the s-p band, which lies well 
below the d states at zero pressure, but well above the lowest d states at P = 300 GPa. We 
have checked that this relative shift is precisely reproduced by the PAW calculations. This 
leaves little doubt that the PAW techniques, on which all our melting calculations are based, 
correctly reproduce this important feature of the pressure dependent electronic structure. 
In Sec. IIVI we will present results on the temperature dependence of the electronic DOS in 
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the solid and the liquid. 

C. Phonon frequencies 

The calculation of phonon frequencies is an important test of DFT approximations, be- 
cause of the detailed comparisons with experimental data that can be made. It is particularly 
important in the context of melting calculations, because, for the harmonic solid, errors in 
phonon frequencies translate directly into free energy errors, which are linked with errors in 
melting temperature. We present here our calculations of the phonon dispersion relations 
of Mo at its experimental volume, using PAW with GGA(PBE) exchange-correlation. 

The technique for calculating the phonon frequencies is the small displacement method, 
as implemented in our phon code^, which we used in earlier work on Fe, Al and Cu. In 
this method, the elements of the force-constant matrix are obtained by displacing atoms 
from the perfect-lattice positions and computing by DFT the forces on all the atoms. For 
a b.c.c. crystal, it suffices to displace a single atom along the (111) direction. Since all 
the calculations employ periodic boundary conditions, the displaced atom is at the centre 
of a periodically repeated supercell. To obtain accurate dispersion relations over the whole 
Brillouin zone, this supercell must be large enough so that the elements of the force- const ant 
matrix have negligible values at its boundaries. In addition, the calculation of the forces must 
be converged with respect to electronic k-point sampling, and to enhance this convergence we 
employ Fermi smearing. This smearing itself incurs errors, which need to be made negligible. 
In summary, the phonon frequencies must be converged with respect to atom displacement, 
supercell size, k-point sampling and Fermi smearing. Rather than insisting that every single 
frequency be converged, it is more convenient to require that the geometric mean frequency 
uj be converged. This quantity is defined by the equation: 



where uj^i is the phonon frequency of branch i at wave vector q, and iV q i is the number of 
branches times total number of q points in the sum. It is useful to work with u, because 
it is directly related to the harmonic free energy of lattice vibrations, which, well above the 
Debye temperature, is equal to 3k^T ln(hu / ksT) per atom. 

Our aim is to have uj converged to better than 1 % with respect to all technical parameters. 




(6) 
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Tests on small supercells show that an atomic displacement of 0.0062 A is small enough to 
ensure that anharmonic errors are well below this tolerance. To test convergence with respect 
to k-point sampling and Fermi smearing width a, we have done extensive tests on a 2 x 2 x 2 
supercell containing 8 atoms (see Table I). These tests show that for a = 0.7 eV, it is easy 
to achieve excellent k-point convergence. But repetition of this at a = 0.5 eV reveals that 
this reduction of a causes u to change by ~ 1 %. However, further reduction of cr to 0.3 eV 
changes u by less than 0.5 %, so that uj appears to be adequately converged with respect to 
k-points and o with a 12 x 12 x 12 k-point grid and a = 0.5 eV. We then seek convergence 
with respect to supercell size using this value of a and a k-point density that is reduced in 
inverse proportion to supercell size. The results indicate that convergence to our required 
tolerance is achieved with the 4x4x4 supercell of 64 atoms. 

Fig. [5] shows a comparison with experimental data of our calculated phonon frequencies 
obtained with the 6x6x6 supercell of 216 atoms, a 4 x 4 x 4 k-point grid, and a = 0.5 eV. 
According to our convergence tests, any discrepancy with experimental frequencies of over 
1 % represents a genuine disagreement. The agreement is actually very satisfactory over 
most of the Brillouin zone, with typical discrepancies being 1 — 2 %. However, there is a 
region around the H point k = (27r/ao)(l, 0, 0), where there are discrepancies of ~ 5 %. This 
same H-point problem has been noted by previous authors who dealt with transition metals 
(Mo and Nb) and used pseudopotential-based methods^ 1 ^. The origin of this sharp dip in 
the phonon dispersion curve of Mo at point H has been related to the nesting of electronic 
states near the Fermi level^ 1 ^ (see figures in Sec. IIII Bl) . so it is likely that by reducing 
the Fermi smearing and increasing the number of k-points the agreement with experiments 
would be improved. In any case, we believe that since the discrepancies are rather localized 
in k-space, they will have only a weak effect on the thermodynamic properties of the system. 

D. Conclusions from the tests 

In summary, our tests show that: (i) neither LDA nor GGA perfectly reproduces the 
experimental T = K pressure-volume curve, but the volume given by GGA deviates by 
only a small and almost constant amount of ~ 1.5 % from the experimental value over the 
pressure range — 300 GPa; (ii) with GGA, the inclusion of 4s states in the valence set 
makes a negligible difference to the P(V) curve; (iii) comparisons of PAW and FP-LAPW 
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confirm the accuracy of PAW for both P(V) and the electronic DOS, and in particular PAW 
accurately reproduces the well-known pressure induced shift of s-p bands relative to d bands; 
(iv) GGA gives rather accurate phonon frequencies over most of the Brillouin zone. This 
evidence provides a firm basis for our calculations on the high-pressure melting of Mo, which 
employ the PAW technique with GGA(PBE) exchange-correlation, and with 4p states but 
not 4s states in the valence set. 



IV. MELTING 



We begin this Section by presenting our ab initio calculations of the melting curve of 
Mo using the reference coexistence technique (see Sec. IIIBj) ; our results for the volume and 
entropy changes on melting as a function of pressure are also reported. In Sec. IIVP4 we 
outline our calculations of the atomic and electronic structures of the b.c.c. solid and the 
liquid. 



A. Calculation of melting curve 

We start by determining the ab initio T m at a pressure close to zero. At this pressure, we 
have experimental values for T m (2883 K)^ and the volumes per atom of coexisting solid 
and liquid (16.34 and 17.04 A/atom) 56 . We use this information to set the temperature and 
volume/atom of the ab initio m.d. simulations that we performed to fit the parameters of 
the reference model. These m.d. simulations employed systems of 125 atoms with T-point 
sampling, and were done at T = 2800 K, V = 16.34 A 3 /atom (solid) and T = 3000 K, V = 
16.34 A 3 /atom (liquid) . The starting configuration of the solid was produced by equilibrating 
the system, which initially was in the perfect crystal configuration, to temperature 2800 K. 
For the liquid, we produced the starting configuration by raising the temperature of the 
perfect crystal to 8000 K (more than twice the experimental T m ) and then rethermalizing 
it again to 3000 K. We checked that the system was in the liquid state by monitoring the 
time-dependent mean-squared displacement. The durations of the simulations were about 
2 ps for the solid and 4 ps for the liquid. A set of 100 solid and liquid configurations from 
these ab initio simulations was then used to fit the reference model by varying the EAM 
parameters to minimize the dimensionless quantity ip and the pressure difference AP (see 
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Sec. Ill Bj) . The resulting EAM parameters are reported in Table II. The very small value 
■if} = 0.078 indicates a good quality of fit, and the AP values of 0.3 GPa and 0.7 GPa for 
the liquid and solid, respectively, were also satisfactory (these AP give contributions of 
only 7 and 3 • 10~ 5 eV/atom to AG). We then performed reference coexistence simulations, 
and found that the solid and liquid remain in stable coexistence over periods of 60 ps at 
P = 11 GPa and T = 3260 K. To illustrate this, we show in Fig. [6] the density profile 
obtained by calculating the number of atoms in slices parallel to the solid-liquid interface. 
The solid is immediately recognisable from the regular oscillations with a repeat distance of 
2.75 A (equal to v^3/2 times the lattice parameter of the bcc lattice, corresponding to the 
distance between nearest neighbours), whereas the density profile is flat in the liquid region. 
Finally, we corrected for the difference between ab initio and reference energy functions, 
obtaining a final ab initio T m = 3205 K at P = 11 GPa. The values of (AUY^/N (where 
{AU)[ s cf = (AU}[ cf - {AU) s rc{ ) and of ((5AU) 2 ) vc{ /2Nk B T of the liquid and solid used to 
make these corrections are reported in Table III. We comment below on the volume and 
entropy of melting. 

We now use the reference model to obtain a first estimate of the melting curve at higher 
pressures. Reference coexistence simulations performed with the EAM parameters from the 
fit at P ~ GPa showed that at P = 92 GPa, the reference T m is 5392 K, the volumes of 
coexisting solid and liquid being 13.30 and 13.53 A 3 /atom. Correcting for the differences 
between ab initio and reference energy functions, we obtain the ab initio melting temperature 
T m = 4867 K at P = 92 GPa. We note that the correction to T m is considerably greater than 
at P ~ 10 GPa, but we believe it is still small enough for the first-order correction scheme 
to remain valid, and this is confirmed by subsequent results (see below). However, when we 
repeated this procedure at P ~ 160 GPa, still using the reference model fitted at P ~ 0, we 
found that the corrections need to go from reference to ab initio T m were even larger, and 
we considered it essential to refit the reference model. Rather than attempting to do this 
refit at P ~ 160 GPa, we returned to P ~ 90 GPa, where our knowledge of the ab initio T m 
is reasonably secure. The refitting at P ~ 90 GPa produced the new parameters reported 
in Table II (100 < P < 200 GPa). This new reference model, when used in coexistence 
simulations at P = 156 GPa, yielded the reference T m = 6510 K, and a corrected ab initio 
value T m = 5969 K. We regarded the size of this correction as acceptably small. In a similar 
way, when we performed calculations at P ~ 270 GPa, using the reference model fitted 
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at 90 GPa, the corrections were unacceptably large, and we performed a refit at 160 GPa. 
This refitted reference model required only rather small corrections when used at 270 and 
380 GPa. 

The reference and ab initio T m as a function of pressure from this full set of calculations 
are reported in Fig. [7J We find that the ab initio T m values can be very well fitted with the 
so-called Simon equation 57 T m = a(l+P/b) c , with a = 2894 K, b = 37.22 GPa and c = 0.433. 
The resulting P = melting temperature of 2894 K is very close to the experimental value 
of 2883 K. Using the Simon equation, we can obtain the melting slope dT m /dP at any 
pressure. At P = 0, we find dT m /dP = 33.7 K GPa -1 , which agrees closely with the 
experimental valued of 33.3 K GPa -1 . Also shown in Fig. [7J is the point on the melting 
curve at P ~ 375 GPa estimated from the shock data of Hixson et a/.-, which is close to 
our melting curve. The diamond anvil cell (DAC) measurements of Errandonea et al— differ 
greatly from our results, since their dT m /dP is essentially zero over most the range from 
to 100 GPa. Previous theoretical melting curves^ for Mo, also shown in the Figure, are 
in general agreement with our results, though there are substantial quantitative differences. 
The comparison of all these experimental and theoretical results raises important issues, 
which will be discussed in Sec. |VJ 

The entropy and volume of fusion of the reference model, denoted by S^ f and V^ e s f , are 
straightforward to calculate. From the reference coexistence simulations, we have (P, T) 
pairs lying on the reference melting curve. We then perform separate single-phase m.d. 
simulations of the solid and liquid reference systems at chosen temperatures, using systems 
of 3375 atoms, adjusting the volumes in each case to give the appropriate P. The difference 
of the resulting volumes give us V^ e s f . At the same time, the simulations give the internal 
energy U, from which we obtain the enthalpies H = U + PV of the two phases. Then the 
difference of enthalpies H\ s ei gives us the entropy difference, since TS\ s e{ = H\ s ei . We find 
that S\ s et is almost constant along the reference melting curve, going from 0.58 /ce/atom at 
P ~ to 0.69 fee/atom at P = 378 GPa. Since the reference system mimics the ab initio 
system rather closely, we assume that SJ^ for the ab initio system is essentially the same as 
S l r l { . To obtain the ab initio volume of fusion V^j, we use the Clausius-Clapeyron relation 
dT m /dP = Vj^l SJ^. Our V ls results for the reference and ab initio systems as a function of 
P are reported in Fig. [HJ We note that in both cases the fractional volume change V ls /V s 
decreases smoothly from ~ 1.5 % at P = to ~ 0.9 % at 400 GPa. 
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B. Atomic and electronic structure of solid and liquid 

We noted in the Introduction that theories of the melting of b.c.c. transition metals 
sometimes assume^ that melting is associated with a significant change of coordination 
number, so that the electronic density of states should also change markedly. In order 
to test these ideas for Mo, we have calculated the radial distribution function g(r) of Mo 
from a series of AIMD simulations of the solid and the liquid performed at P ~ 216 GPa. 
The simulations were all done with 125 atoms, and had a typical duration of 2 ps after 
equilibration. Fig. [9] reports our calculated g(r) at T = 2000, 4000 and 6000 K (solid state), 
and at T = 7500 K (liquid state); we recall (Fig. [7]) that our calculated T m at this pressure 
is 6641 K. At 2000 K, the shells containing 8 first neighbours and 6 second neighbours at 
distances of r = 2.44 and 2.81 A are clearly separated, and g(r) goes to zero at r ~ 3.4 A, 
between the second and third shells. However, at T = 4000 K, the first and second shells 
have already merged, though the shoulder due to the second shell is clearly visible. The 
g(r) of the solid at T = 6000 K presents similar trends to those shown at 4000 K, the 
main difference being that the valley between the second and third atomic shells now clearly 
moves upwards from zero. The change in g(r) in going from the solid at T = 6000 K to 
the liquid at 7500 K is substantial for r > 3 A, with the well defined peaks due to third 
and higher neighbours in the solid becoming heavily broadened. However, the peak closest 
to the origin does not suffer a large change. We define the coordination number N c in the 
usual way as N c = 4irp g{r)r 2 dr, with p the bulk number density and r c the distance 
of the first minimum. The r c values at the four temperatures are 3.41, 3.41, 3.41 and 
3.31 A, and the resulting N c values are 14, 14, 14 and 13.35. Apart from the expected 
increase in disorder, the main change on going from solid to liquid is thus a slight decrease 
in coordination number. We comment further on this in Sec. PVT . 

Turning now to the electronic density of states (DOS), we present first our AIMD results 
for P in the range 50 — 70 GPa at a series of temperatures, the simulations being performed 
on a system of 125 atoms with T-point sampling. The typical duration of these simulations 
was 2 ps after equilibration and the DOS were calculated by averaging over 150 different 
configurations. We report in Fig. [10] the calculated DOS at the thermodynamic states given 
by (P, T) = (48, 0), (50, 2000), (51, 3300) for the solid, and (P, T) = (72, 5000) for the liquid 
(units of GPa and K). (Our ab initio T m for 50 < P < 70 GPa are in the range 4185—4577 K.) 
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We have checked in each case that the system is in the solid or liquid states by looking at 
the mean-squared displacement and structure factor. We note the progressive broadening 
of the DOS peaks with increasing thermal disorder in the solid, an effect which continues 
further in the liquid. The Fermi-level value of the DOS increases slightly on melting. As far 
as occupied states are concerned, melting appears to cause a slight redistribution of d-states 
from lower in the band to the region of the Fermi level. 

In the right-hand panel of Fig. [TOj, we compare our AIMD results for the electronic DOS 
at the solid state-point (P, T) = (285, 7000) and the liquid state-point (P, T) = (300, 8250) 
(units of GPa and K), which are just below and just above our calculated melting curve. We 
note the rather minor changes caused by melting. Interestingly, the Fermi-level value of the 
DOS is almost identical in the two phases at this pressure. The relationship of these results 
with earlier work on the electronic structure of liquid transition metals will be discussed in 
the following Section. 



V. DISCUSSION AND CONCLUSIONS 



At the start of this paper, we emphasized the large discrepancies between melting curves 
of transition metals derived from static compression and shock measurements, and we men- 
tioned that previous DFT work on Mo supports the shock measurements. The present 
work fully confirms that the melting curve predicted by DFT in the PBE approximation for 
exchange-correlation energy lies far above the static compression measurements, but at high 
pressures is consistent with the shock data. This confirmation is important, because of defi- 
ciencies or uncertainties in previous DFT work. The reliability of the present calculations is 
supported by our close agreement with the experimental P = values of both the melting 
temperature T m and the melting slope dT m /dP. Our melting curve is below the theoretical 
curve of Moriarty 6 by ~ 600 K at P = 0, and this difference increases with increasing P. 
However, this is not surprising, since the generalized pseudopotential model that he used is 
known to disagree with experimental phonon frequencies, and because he included thermal 
electronic excitations only approximately. In the present work, we have taken pains to verify 
the accuracy of the phonon frequencies given by our methods, and thermal electronic exci- 
tations are fully included within our DFT framework. Perhaps more surprising is that our 
melting curve agrees closely with that obtained by Belonoshko et air- using direct DFT m.d. 
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simulation. This is unexpected, since they believed that their melting curve suffered from a 
substantial superheating error of ~ 20 %. The close agreement suggests that they may have 
been unduly pessimistic, and this point deserves further investigation. We included in Fig. [7] 
the results of Belonoshko et al— and of Verma et al— obtained by the dislocation-mediated 
theory of melting^^. It is not clear to us whether one can expect a theory of melting based 
exclusively on the properties of the solid to be fully reliable. One of the problems with this 
approach is that the predicted melting curves rely on thermodynamic data that may not be 
reliably known. The rather large differences between the two melting curves based on the 
dislocation theory may be indicative of the limited reliability of this approach. 

The change of volume on melting of ~ 1 % given by our calculations is small, but still 
much greater than the volume change implied by the static compression values^^i of 
dTja/dP. Arguments in favour of a very low volume change based on a significant increase 
of coordination number on going from b.c.c. solid to melt appear to be incorrect, according 
to our DFT m.d. calculations of the radial distribution function g(r). We find only rather 
minor differences between g(r) for high-T solid and melt. In particular, there is actually 
a slight decrease in coordination number from 14 to ~ 13 on melting, so that the liquid 
is slightly less close packed than the solid. We comment that ideas based on hard-sphere 
packing are likely to be misleading, since the repulsive interactions between Mo atoms at 
high T are rather soft (see below). 

We mentioned in the Introduction the recent theory of Ross et al.—, according to which 
a very low melting slope is expected for b.c.c. transition metals. The theory invokes the 
well known transfer of electrons from s-p to d states with increasing compression, and the 
fact that this transfer depends on crystal structure. In applying this theory to the melting 
of Mo, the authors estimated the effective number of d electrons n& by treating the high- 
temperature solid as a perfect b.c.c. crystal and the liquid as a perfect f.c.c. crystal. They 
also assumed that a change of on melting will be associated with a change of d-band 
width. They found that the changes of electronic structure stabilize the liquid relative to 
the solid, and yield a major reduction of T m . In considering this theory in the light of the 
results we have presented, it is important to appreciate that our calculations are all based on 
an accurate implementation of DFT. As in all simulations using Born-Oppenheimer DFT 
m.d., the VASP code recalculates the entire self-consistent electronic structure at every 
time step of the time evolution. As described in Sec. IIII Bl we have gone to considerable 
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lengths to show that the PAW implementation of DFT used in our m.d. yields results for 
the electronic DOS which are almost indistinguishable from those given by the FP-LAPW 
technique, which is one of the most accurate available. This means that all the effects that 
enter the theory of Ross et al.—, including s-p to d electron transfer and changes of d-band 
width, are fully included in our simulations. Nevertheless, we do not obtain the very low 
melting slope that they predict. The reason for this is presumably that their treatment 
of the high-T solid and the liquid as perfect crystals is incorrect. As we have seen, their 
assumption of a large structural change on melting also appears to be questionable. This 
point is reinforced by our finding that the electronic DOS changes only slightly on melting, 
especially at high P. 

We end this discussion by commenting on the embedded-atom model (EAM) used as a 
reference system in determining the melting curve. It is an important finding of this work 
that the EAM is able to mimic accurately the DFT total-energy function of solid and liquid 
close to coexistence. This does not mean, however, that we accept the melting curve of 
the reference model as the true melting curve. This could be dangerous, because we might 
then miss d-band electronic effects that were not explicitly included in the model. However, 
in our procedure, any such effects are automatically picked up in the corrections that we 
apply, since these explicitly account for free energy differences between the reference and 
DFT systems. We note in passing that the fitting of our reference model yields parameters 
that resemble those we found in our earlier work on the melting of Fe^. In particular, the 
inverse-power repulsive potential in our present EAM model has an exponent n close to 6 
at low P, decreasing to ~ 5 at high P. For comparison, the fitted EAM in our Fe work 
had n = 5.9, which is very similar. We are currently investigating the systematic behaviour 
of EAM parameters in solid and liquid transition metals at high P and T, and we hope to 
report on this elsewhere. 

The main conclusions from this work are as follows. Our DFT calculations of the melting 
curve of Mo up to 400 GPa fully confirm earlier work showing that the DFT melting curve 
is consistent with shock data, but is far above the melting curve given by static compression 
experiments. Our calculations indicate that at high P there are only minor changes of 
both atomic and electronic structure on going from the high-temperature b.c.c. solid to the 
melt. Suggested mechanisms for an anomalously low melting slope dT m /dP of Mo based on 
transfer of electrons from s-p states to d states appear to be incompatible with the present 
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DFT calculations. This tends to confirm earlier suggestions 9 that the transition identified as 
melting in high-P static compression experiments may not be true thermodynamic melting. 
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N 


supercell 


k-grid 


a (eV) 


57 (10 12 s" 1 ) 


8 


2x2x2 


4x4x4 


0.7 


20.734 






8x8x8 


0.7 


21.385 






12 x 12 x 12 


0.7 


21.204 






16 x 16 x 16 


0.7 


21.206 


8 


2x2x2 


8x8x8 


0.5 


21.498 






12 x 12 x 12 


0.5 


21.045 






16 x 16 x 16 


0.5 


21.037 






24 x 24 x 24 


0.5 


21.045 


8 


2x2x2 


8x8x8 


0.3 


21.703 






12 x 12 x 12 


0.3 


20.921 






16 x 16 x 16 


0.3 


20.811 






24 x 24 x 24 


0.3 


20.951 






32 x 32 x 32 


0.3 


20.941 


64 


4x4x4 


6x6x6 


0.5 


21.699 


216 


6x6x6 


4x4x4 


0.5 


21.727 



TABLE I: Convergence of mean phonon frequency uj (see Eq. [6]) with supercell size, k-grid and 
Fermi broadening a. 
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P (GPa) e (eV) a (A) n m C 



- 


- 100 


0.323 


3.579 


5.93 


3.72 


12.66 


100 


- 200 


0.169 


4.985 


4.96 


3.88 


28.08 


200 


-400 


0.144 


4.760 


5.07 


3.78 


26.90 



TABLE II: Parameters of the EAM potential deduced for Mo and used for the coexistence simu- 
lations. Values are obtained by fitting to ab initio simulations on solid and liquid. 
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T™ { (K) 


(AU) l JLJN (eV/atom) 

\ / rCI f \ / / 


y((SAuf) rci 


■MT (eV/atom) 


T^ 1 (K) 






Solid 


Liquid 




3260 


0.009(2) 


0.038(2) 


0.032(2) 


3205 


5392 


0.027(1) 


0.024(2) 


0.028(1) 


4867 


6510 


0.038(1) 


0.035(3) 


0.030(2) 


5969 


7324 


0.002(1) 


0.023(2) 


0.015(2) 


7281 


8618 


0.013(1) 


0.018(2) 


0.032(2) 


8154 



TABLE III: Difference {AU) l r s ef = (A?7)[ ef — {AU)^ el between the liquid and solid thermal averages 
of the difference AU = Uai — U rc f of ab initio and reference energies, and thermal averages in 
solid and liquid ((SAU) 2 ) rc f of the squared fluctuations of SAU = AU — (AU) rc {, with averages 
evaluated in the reference system and normalized by dividing by the number of atoms N. Melting 
temperatures for the reference and ab initio systems are also reported. 
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Figure Caption List 



FIG. [TJ Comparison of LDA and GGA pressure P as function of volume V for b.c.c. 
Mo from the PAW method with different exchange-correlation functionals and valence sets. 
Long-dashed and solid lines (practically coincident) show GGA results with and without 4s 
states in the valence set. Short-dashed and dotted lines show LDA results with and without 
4s states in the valence set. Dots show experimental results 5 . 

FIG. [2J Comparison between PAW and FP-LAPW results for the GGA(PBE) and 
LDA(CA) approximations for E xc . Solid and dashed curves show GGA(PBE) and LDA(CA) 
FP-LAPW results, respectively; short-dashed and dotted curves show GGA(PBE) and 
LDA(CA) PAW calculations, respectively. Solid dots show experimental data 5 . 

FIG. EJ Density of electronic states obtained with the PAW (dashed line) and FP-LAPW 
(solid line) at zero temperature and 0, 150 and 300 GPa. Fermi energies are shifted to zero 
(dotted line). 

FIG. SJ FP-LAPW calculation of the energy bands of Mo at and 300 GPa (left and 
right panels respectively). The s valence band (energy between -8 and -7 eV at the V point) 
rises in energy more quickly than d valence bands with increasing pressure. 

FIG. [5j Comparison of calculated (curves) and experimental (solid squares) phonon 
dispersion relations of Mo at zero pressure. Experimental data are from Ref.— . 

FIG. [6j Density profile in simulation of coexisting solid and liquid Mo at P = 11 GPa, 
T = 3260 K after 60 ps. The simulation is performed with the embedded-atom reference 
model on a system of 6750 atoms. 

FIG. [7J Calculated ah initio melting curve (filled circles and solid line) of this work com- 
pared with previous results: generalized pseudopotential calculations of Moriarty 6 (dotted 
line), dislocation-mediated models of Belonoshko et al. 9 (long-dashed line) and Verma et 
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al— (dashed-dotted line); experimental shock- wave 4 and DAC^ measurements are shown 
with empty squares and triangles, respectively. Filled and inverted-empty triangles show 
solid and liquid ab initio molecular dynamics calculations of Belonoshko et al. 9 , respectively. 
Empty circles show results of this work obtained with the EAM model without free-energy 
corrections. 

FIG. [HJ Ab initio fractional volume change on melting of Mo as a function of pressure. 
Solid and dashed curves: present work, with and without free-energy correction, respectively. 

FIG. IHJ Calculated radial distribution function of Mo for: solid at P = 212 GPa and 
T = 2000 K (solid line), solid at P = 216 GPa and T = 4000 K (dotted line), solid 
at P = 230 GPa and T = 6000 K (long-dashed line) and liquid at P = 224 GPa and 
T = 7500 K (short-dashed line). 

FIG. HOI Density of valence electronic states of Mo at finite temperature and on melting. 
Left: solid at P = 48 GPa and T = K (dotted line), solid at P = 50 GPa and T = 2000 K 
(short-dashed line), solid at P = 51 GPa and T = 3300 K (solid line) and liquid at P = 72 
GPa and T = 5000 K (long-dashed line). Right: solid at P = 285 GPa and T = 7000 K 
(solid line) and liquid at P = 300 GPa and T = 8250 K (dashed line). Fermi energy levels 
are shifted to zero. 
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FIG. 2: 
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FIG. 4: 
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FIG. 9: 
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